Journal of Glaciology, Vol. 59, No. 216, 2013 doi:10.3189/2013JoG12J147 


613 


Antarctica, Greenland and Gulf of Alaska land-ice evolution from 
an iterated GRACE global mascon solution 

Scott B. LUTHCKE , 1 T.J. SAB AKA , 1 B.D. LOOMIS , 2 A.A. ARENDT , 3 J.J. McCARTHY , 2 

J. CAMP 4 

1 Planetary Geodynamics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA 

E-mail: scott.b.iuthcke@nasa.gov 
2 Science Division, SGT Inc., Greenbelt, MD, USA 
3 University of Alaska, Fairbanks, AK, USA 

4 Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA 

ABSTRACT. We have determined the ice mass evolution of the Antarctic and Greenland ice sheets (AIS 
and GIS) and Gulf of Alaska (GOA) glaciers from a new GRACE global solution of equal-area surface 
mass concentration parcels (mascons) in equivalent height of water. The mascons were estimated 
directly from the reduction of the inter-satellite K-band range-rate (KBRR) observations, taking into 
account the full noise covariance, and formally iterating the solution. The new solution increases signal 
recovery while reducing the GRACE KBRR observation residuals. The mascons were estimated with 
1 0 day and 1 arcdeg equal-area sampling, applying anisotropic constraints. An ensemble empirical mode 
decomposition adaptive filter was applied to the mascon time series to compute annual mass balances. 

The details and causes of the spatial and temporal variability of the land-ice regions studied are 
discussed. The estimated mass trend over the total GIS, AIS and GOA glaciers for the time period 
1 December 2003 to 1 December 2010 is -380 ±31 Gta -1 , equivalent to -1.05 ±0.09 mm a -1 sea-level 
rise. Over the same time period we estimate the mass acceleration to be -41 ± 27 Gta -2 , equivalent to a 
-0.11 ±0.08 mm a -2 sea-level acceleration. The trends and accelerations are dependent on significant 
seasonal and annual balance anomalies. 


INTRODUCTION 

The land-ice component of the cryosphere has important 
connections to the global climate system. Changes in land- 
ice cover impact the Earth's radiation budget through 
variations in albedo, and significantly impact the hydro- 
logical cycle as changes in water storage. The Earth's large 
areas of land ice, including the ice sheets, glaciers and ice 
caps, have the potential to significantly impact eustatic sea 
level, as they hold ~77% of the world's fresh water, which 
amounts to ~65 m of global sea-level equivalent (Barry and 
Gan, 2011). Changes in freshwater input from calving and 
melting Arctic land and sea ice can increase the haline 
forcing in the North Atlantic and impact global ocean 
circulation and climate patterns (Hakkinen and Rhines, 
2004). Variations in land-ice volume are closely correlated 
to global temperatures, as paleoclimate records show all the 
major Greenland ice sheet (GIS) changes have been 
correlated with temperature changes (Alley and others, 
201 0). It is important that we accurately quantify the current 
spatial and temporal mass changes of the Earth's land ice to 
improve our understanding, modeling and prediction of 
its response and contribution to climate change. Since its 
launch in March 2002, the Gravity Recovery and Climate 
Experiment (GRACE) mission has acquired ultra-precise 
inter-satellite K-band range and range-rate (KBRR) measure- 
ments between two co-orbiting satellites in a 450km 
altitude polar orbit, ~220km apart (Tapley and others, 
2004). These data have vastly improved our knowledge of 
the Earth's time-variable gravity field. In particular the 
GRACE mission has been instrumental in quantifying recent 
land-ice mass evolution (e.g. Luthcke and others, 2006a, 
2008; Velicogna, 2009; Chen and others, 2011; Schrama 


and Wouters, 2011; Jacob and others, 2012; King and 
others, 2012). 

There are numerous challenges in applying GRACE to the 
study of glacier and ice-sheet mass balances. GRACE senses 
all sources of mass variability at the Earth's surface, requiring 
removal of non-glacier changes through incorporation of 
independent datasets and models. Mass signals associated 
with glacial isostatic adjustments (GIA) and terrestrial water 
storage (TWS) are poorly known in many regions and often 
dominate GRACE error budgets. Owing to the GRACE 
mission's limited longitudinal sampling at monthly time 
intervals, filtering techniques are often applied to mitigate 
striping errors and to focus signal into a region of interest, 
often resulting in signal attenuation. When mass is 
constrained to a particular spatial extent, signal can either 
leak into or out of that domain and cause additional error. 
Even when the choice of filtering and models used to correct 
atmosphere and ocean mass variations is fixed, differences 
in GRACE solutions can occur, due to different methods 
used in processing the fundamental GRACE observations 
(Pritchard and others, 2010). 

Here we present a new study to quantify the ice mass 
evolution of the Antarctic and Greenland ice sheets (AIS and 
GIS) and Gulf of Alaska (GOA) glaciers from a GRACE 
global mascon solution. The mascon technique is a time/ 
space domain approach for describing time-variable gravity 
signals, in which time is discretized in contiguous, non- 
overlapping finite intervals and, within a given interval, 
space is divided into gravity contributions due to loading of 
mass concentrations of equivalent water (mascons) on the 
Earth surface. This study is unique and differs in many 
aspects from previous regional and global mascon solutions, 
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including those of Luthcke and others (2006a, 2008), 
Rowlands and others (2010) and Sabaka and others 
(2010). The mascons are estimated globally with 10 day 
temporal and 1 arcdeg equal-area spatial sampling. Ice- 
sheet mascons are restricted to the grounded-ice areas 
contributing to eustatic sea-level change. We apply tem- 
poral and spatial constraints anisotropically, through seven 
distinct constraint regions representing the ice sheets (high- 
interior and low-elevation margins), GOA glaciers, land and 
ocean. The constraint equation is an exponential function of 
correlation time and distance applied to all possible pairs of 
mascons within a constraint region (Sabaka and others, 
2010). A detailed error analysis is provided, including 
consideration of leakage, and errors in modeling mass 
variations of the atmosphere and ocean during the GRACE 
data reduction (forward modeling). In addition, a controlled 
comparison is performed in which we use the same set of 
forward modeling and GRACE data reduction, but compute 
the mass time series using an averaging kernel filter applied 
to monthly spherical harmonic fields, as is common 
practice in the literature (e.g. Velicogna and Wahr, 2005). 
The mascon solution is iterated until convergence, and the 
impact of the iterated solution on the reduction of the 
GRACE KBRR residuals is quantified. Furthermore, this study 
differs from the mascon solution presented by Jacob and 
others (2012), in that we estimate the global mascons 
directly from the GRACE KBRR data, taking into account the 
full noise covariance, and iterate the solution while 
quantifying the minimization of the GRACE KBRR obser- 
vation residuals. We present the details of the global mascon 
solution, an error analysis and a discussion of results focused 
on recent GIS, AIS and GOA glaciers land-ice evolution. 
This paper also serves to document our mascon solution 
used in the international Ice Mass Balance Inter-comparison 
Exercise (IMBIE; Shepherd and others, 2012). 


MASCON FORMULATION 

The GRACE processing centers provide time-variable gravity 
observations (Level-2 products) in the form of monthly sets 
of Stokes coefficients or spherical harmonic fields. These 
monthly fields are estimated directly from GRACE tracking 
data, typically without the use of constraints or a priori 
information. The use of these fields for scientific applica- 
tions, such as measuring the mass evolution of the GIS or 
hydrological basins, requires filtering due to north-south 
striping artifacts, which arise from the limited longitudinal 
sampling of the GRACE orbits within monthly time periods 
(Swenson and Wahr, 2006; Klees and others, 2008). Many 
filtering methods do not take into account the noise 
covariance of the estimated Stokes coefficients, and there- 
fore are not optimal filters (Klees and others, 2008). Loss of 
signal amplitude and limited spatial and temporal reso- 
lution, as well as signal leakage in and out of regions of 
interest, are particular problems when applying filtering 
methods. 

Instead of applying a filtering technique to the GRACE 
Level-2 monthly Stokes coefficients, we used a mascon 
technique, applying geolocatable anisotropic constraints to 
estimate the global mass change directly from GRACE KBRR 
data, which accounts for the full Stokes noise covariance. 

The formulation for mascon parameters is derived from 
the fact that a change in the gravitational potential caused by 
adding a small uniform layer of mass over a region at an 


epoch t can be represented as a set of (differential) potential 
coefficients which can be added to the mean field. The delta 
coefficients can be computed as (Chao and others, 1987) 
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where k\ is the loading Love number of degree /, to account 
for the Earth's elastic yielding which, in general, counteracts 
the additional surface density (Farrell, 1972); A, cp and R are 
geographic latitude, geographic longitude and the Earth's 
mean semi-major axis, respectively; M is the Earth's mass; 
H(t) is equivalent water height in centimeters of the mass of 
the layer over a unit of surface area at the epoch t; and fi is 
the solid angle surface area of the mascon region where H(t) 
is applied, dfl = cosAdAd^. 

The water height, /-/(f) (cm), can be equated to a surface 
mass density, a(t) (kgm -2 ), by noting that water has a 
volume density of 1000 kgm -3 , which gives <r(f) = 
10 x H(t). The estimated mascon parameter, /-/,■( f), for each 
mascon region, j, is a scale factor on the set of differential 
Stokes coefficients for that mascon region, giving a surface 
mass change (cmw.e.). In matrix notation, an assemblage 
over j of the above equations provides a set of partial 
derivatives, L, of the differential Stokes coefficients, As, 
(complete to degree and order 60) with respect to equivalent 
water height, Ah, such that 

As = LAh. (3) 

The mascon parameters are part of nonlinear functions of 
the GRACE KBRR measurements and are therefore estimated 
using the Gauss-Newton (GN) method (Seber and Wild, 
1989) of nonlinear least-squares. This method seeks a series 
of differential corrections to the parameters via a lineariza- 
tion of the observation equations at the current state. This 
means that the changes in KBRR observations at the current 
state due to changes in mascon parameters are needed for 
the estimation. The partial derivatives of the GRACE KBRR 
observations with respect to the mascon parameters are 
computed by the chain rule, using the partial derivatives of 
the KBRR observations with respect to standard Stokes 
coefficients and the differential Stokes coefficients com- 
puted from the equation above for the mascon region of 
interest: 


dQi ' dOjdACijt) dO,dAS j lm (t) 

tr^o dC im dH j(t) dS lm dHj(t) ' 

where dOi/dHj(t) is the partial derivative of KBRR obser- 
vation / with respect to the j mascon parameter, Hy(f), at 
time f; dO\/dCi m and <90,/3S/ m are the partial derivatives of 
the KBRR observations with respect to the geopotential 
Stokes coefficients (computed as part of the KBRR reduction 

and Level-1 data processing); and d&.C\ m {t)/dHj(t) and 

9AS ; /m (f) / dHj(t) are the partial derivatives of the delta 
Stokes coefficients with respect to the mascon parameter in 
region j at time t. 

In matrix notation, assemblages over / and j of the above 
equation provide a set of partial derivatives of the KBRR 
observations, o, with respect to the mascon parameters, h, 
by chaining L with partial derivatives, A, of the KBRR 
observations with respect to the Stokes coefficients. 
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Following Sabaka and others (2010), at step k of the GN 
estimation let the GRACE KBRR residuals, r, be defined as 
the difference between the observations, o, and a model 
prediction, a(x,), which is a function of a current state of a 
general parameter set, x,. In addition to the current mascon 
parameter state, this parameter set represents the effects of 
orbital arc parameters (initial baseline state orbit and 
accelerometer calibration parameters) and the mean field 
and forward models of other geophysical processes of mass 
redistribution (e.g. atmosphere and ocean tides). The details 
of the orbit arc parameters and forward models are discussed 
below, in the data processing section. The correction, Ah*., 

to the current state of mascon parameters, h,, is then defined 
as the unique minimizer of the quadratic cost function 

J( Ah,) =(r - ALAh,) T W(r - ALAh,) 

+ A*(h, + Ah,^ P/,/i (h k + Ah,^ , 


where the first term is the weighted misfit and the second 
term is a measure of mascon complexity, which in this case 
is a smoothing function that tries to drive the mascon state to 
zero. The data weight matrix, W, accounts for both 
measurement noise and the effects of the orbital arc 
parameters (back substituted), and P,, is a mascon 
regularization matrix, whose influence is controlled by the 

damping parameter, fi. The updated mascon state, h, +1 , is 
given by 


Ah, 


GN iteration k< 


(L t A t WAL + M P hh )~ 
(L T A T Wr - /xPji,) 


. h, +1 = h, + Ah, 


( 6 ) 


assuming a step size of unity. 

The cost function of Eqn (5) suggests a major philosoph- 
ical difference between this study and past studies, in that 
multiple iterations of GN are actually carried out. For 
instance, Rowlands and others (201 0) and Sabaka and others 

(2010) carry out a single step of GN, in which hi is zero in 
Eqn (5). In this study, the step k update, x, +1 , which reflects 

h, +1 , is fed back into the KBRR processing in order to carry 
out the next iteration. In this sense, the minimum of 7(h) is 
sought instead of the local quadratic approximation to it at 

h,, i.e. 7(Ah,). It will be shown below that this leads to a 
significant improvement in the recovered mascon solutions, 
particularly in the ability to restore more signal strength. 

There are several ways to rewrite the GN iteration k that 
illuminate different properties of the estimator. For instance, 
the filtering nature of the estimator and its error sources can 
be seen by first expressing Eqn (6) as 


GN iteration k< , k t\ ,_i T T / - \ 

\ (L t A t WAL + uPm) L T A T w(r + ALh, J 

( 7 ) 

and then assuming the residuals can be expressed as a first- 
order linear perturbation, Ah,, about a true mascon state, 
h, = h,, with additive noise, e, such that r = ALAh, + e. 
This gives 


h,+i = Rh, + Ke, (8) 

where K = (L T A T WAL + /jP,,) -1 L t A t W and R = KAL. The 
estimate, h, +1 , differs from the true state, h, +1 = h, + Ah,, 


due to two effects: (1) the degradation of resolving power in 
the filter, due to the presence of regularization, which is 
manifested as a deviation of the resolution matrix, R, from 
an identity matrix in the first term of Eqn (8), and (2) the 
presence of additive noise from the second term in Eqn (8). 
The effects of the resolution will be discussed in detail 
below. 

Another way to rewrite the GN iteration k is 


GN iteration 


s, = Lh, 

As, = (A T WA) -1 A T r 

s ,+1 — s , + As, 


. h, +1 = P^ L t [ m (A t WA) _ 1 +LP^ L t ] s,+i . 


( 9 ) 


This is now in the form of an anisotropic, non-symmetric 
(ANS) filter (Kusche, 2007) that acts upon an unfiltered 
Stokes state, s, +1 , to produce a filtered Stokes state, 5, +1 , 

~s, + i =P- 1 (N" 1 +P- I )“ 1 s, +1 , (10) 

followed by a least-square co-location (LSC) prediction 
(Moritz, 1980) of the mascon state, h, +1 , from the filtered 
Stokes state 


h , +1 = PjP ss S,+i, (11) 

where N -1 and P" 1 are the noise covariance and signal 
auto-covariance matrices of the Stokes coefficients, respect- 
ively, and Pjjj 1 is the signal cross-covariance matrix between 
the mascons and Stokes coefficients. The ability to process 
the KBRR tracking data directly allows for access to the 
actual noise covariance matrix of the Stokes coefficients, 
and thus makes this an example of an optimal ANS filter 
(Kusche, 2007). Note that the Stokes signal covariances are 
related to the mascon signal auto-covariance through simple 
propagation, such that P'^LPj^L 7 and Pjjj 1 = P kk L T . 
Currently, the mascon signal auto-covariance matrix reflects 
a first-difference operation between all distinct pairs of 
mascons in space and time, weighted as an exponential 
function of temporal and spatial distance between them, and 
a conservation of mass constraint that keeps the global 
mascon sum constant over time. Because of the regional 
aspect of the spatial weighting, the gravity signal is treated as 
a non-stationary, anisotropic process based upon geolocat- 
able physical properties. Details of the constraint regions, 
mascon signal auto-covariance construction and selection of 
the damping parameter are given in the following section. 


GRACE DATA PROCESSING, AND MASCON 

SOLUTION DETAILS 

Data processing and forward models 

We estimated global surface mass anomalies, or mascons, 
directly from the reduction of the GRACE Level-1 B data, 
which includes the inter-satellite K-band range-rate (KBRR) 
data, as well as orbit position, attitude and accelerometer 
data for each GRACE satellite. We employed a baseline state 
parameterization (Rowlands and others, 2002) and accel- 
erometer calibration (Luthcke and others, 2006b), which 
enables the estimation of surface mass variability from 
GRACE KBRR data alone. The GPS data are used solely to 
establish an accurate orbital reference and for calibrating 
accelerometers (Luthcke and others, 2006b). The GRACE 
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Table 1 . Summary of forward modeling applied for each global 
mascon solution 


Mascon solution 

Forward modeling used in computation of KBRR 
residuals number and partial derivatives for 
indicated mascon solution 

v08 

CCM02C (150 x 150) 

Pole tide (IERS2003) 

Solid Earth tides (IERS2003) 

GOT4.7 ocean tides (Ray, 1999; Ray and Ponte, 
2003) 

OMCT ocean model from AOD1 B 
ECMWF 3 hour operational pressure grids 
(N max = 90) 

Post-solution geocenter correction (Swenson and 
others, 2008) 

v09 

Same as v08 

CLDAS/Noah 0.25°, 3 hour data (N max = 90) Ice 
'masked out' 

ICE-5C (VM2) (model C 2 , and S 2 , from GIA/G RACE 
solution)* 

vIO 

Same as v09 plus v09 global mascon solution 

vll 

Same as v09 plus vIO global mascon solution 

vl 2 

Same as v09 plus vll global mascon solution 


•This CIA model is restored in the mascon solutions, and more recent GIA 
models noted in the text are used in the computation of the mascon results. 


KBRR and Level-IB attitude and accelerometer data (in- 
cluding accelerometer calibration) are processed in daily 
arcs as outlined by Luthcke and others (2006b). Forward 
models for both static and time-varying gravity are applied in 
the KBRR data reduction daily arcs, in order to isolate land- 
ice signal. Static mean gravity is modeled using the 
complete GGM02C GRACE gravity model through degree 
and order (d/o) 150 (Tapley and others, 2004). Ocean tides 
are modeled with the GOT4.7 tide model (Ray, 1999; Ray 
and Ponte, 2003). Atmospheric mass variations are modeled 
to d/o 90 at 3 hour intervals using European Centre for 
Medium-Range Weather Forecasts (ECMWF) operational 
pressure grids. Ocean mass variations are modeled to d/o 90 
from the baroclinic Ocean Model for Circulation and Tides 
(OMCT). Terrestrial water storage (TWS) variations are 
modeled using the Global Land Data Assimilation System 
(GLDAS)/Noah Land Surface Model (Noah) 0.25°, 3 hour 
data, represented as potential coefficients to d/o 90 (Ek and 
others, 2003; Rodell and others, 2004). The TWS contri- 
bution was set to zero for the major mountain glacier areas 
of the globe using a 0.25° mask (Raup and others, 2000) and 
also set to zero for the GIS and AIS. Glacial isostatic 
adjustments are modeled based on ICE-5G (VM2) (Peltier, 
2004) and are represented as potential coefficients to d/o 90. 
Table 1 provides a summary of the modeling applied in the 
reduction of the KBRR residuals for each of the mascon 
solutions noted. 

GIA and geocenter 

We modified the mascon solution to restore the forward- 
modeled global GIA, while more recent GIA models are 
then used in the computation of the mascon solution results 
presented. For both the GIS and GOA, a GIA model based 
on the ICE-5G deglaciation history and an incompressible 
two-layer approximation to the VM2 viscosity profile are 
used to correct the GRACE mascon solution (Paulson and 


others, 2007, computed and provided by A. Geruo, 2011). 
For the AIS, the IJ05 _R2 regional GIA model is used (Ivins 
and others, in press). We correct our mascon solutions for 
the viscous component of post-Little Ice Age (post-LIA) GIA 
following the collapse of the Glacier Bay Icefield. As in 
Luthcke and others (2008), we apply a regional post-LIA GIA 
model that is rigidly constrained by ~100 precise GPS and 
relative sea-level (RSL) observations of uplift (Larsen and 
others, 2005). We further correct our mascon solutions to 
reflect surface mass variability in the Earth's center-of-figure 
frame, rather than the center-of-mass frame in which the 
solutions are performed. We apply the geocenter correction 
used in the IMBIE study (Shepherd and others, 2012) derived 
from the degree 1 Stokes coefficients determined from 
Swenson and others (2008). The geocenter correction 
accounts for ~1% of the GIS and GOA ice trends, and 
<3% for the West AIS and AIS peninsula ice trends. The 
largest impact of the geocenter correction is for East AIS, at 
1 7% of the trend. 

Global mascon definitions 

We estimate a global set of 1 x 1 arcdeg equal-area 
(~1 2 390 km 2 ) mascons at a 10 day temporal sampling. 
The choice of the spatial sampling is a compromise between 
approximating regional shapes and boundaries (e.g. land, 
ocean, land ice) and the total number of parameters that can 
be reasonably iterated for a global solution. There are a total 
of 41 168 mascons estimated every 10 days, representing a 
significantly smaller spatial sampling than the fundamental 
spatial resolution of GRACE (Luthcke and others, 2006b). 
Furthermore, the number of mascon parameters is signifi- 
cantly larger than the number of d/o 60 Stokes coefficients 
forming our basis functions (3718). Therefore, as discussed 
in detail above, we employ a regularization in the form of 
the mascon signal auto-covariance matrix (Sabaka and 
others, 2010). This matrix is constructed from constraint 
equations which require all mascon differences to be close 
to zero in a statistical sense. The constraints are applied in 
groups that pertain to geographical regions, so if two 
mascons reside within the same region, the constraint 
weight is a simple exponential function of the distance 
(characteristic length) and time (characteristic time) between 
the two mascons. If two mascons reside in different regions, 
the weight is zero and the constraint vanishes. The details of 
the constraints and the construction of the regularization 
matrix, P^,, are given by Sabaka and others (2010, their 
section 2.2). The important implementation parameters are 
the characteristic length and time used in the exponential 
taper weighting of the mascon first-difference constraint 
equations, and the damping parameter, which controls the 
influence of the regularization relative to the data misfit. 

Application of regional constraints 

The constraints are then applied anisotropically using the 
following seven constraint regions: (1) GIS margins defined 
as <2000 m elevation, (2) GIS >2000 m elevation, (3) AIS 
<2000 m elevation, (4) AIS >2000 m elevation, (5) GOA 
glacier area, (6) land (including all other ice regions) and 
(7) ocean and large bodies of water (including ice shelves). 
In defining the mascons for land and large bodies of water, a 
minimum enclosed region area of 80 000 km 2 was used. 
Thus, if a land area enclosed within the ocean region is 
<80 000 km 2 , then those land mascons were labeled as 
ocean (similarly for water regions enclosed within land). 
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Fig. 1 . Land-ice mascon configurations and drainage systems, (a) The GIS. The numeric labels describe the primary drainage systems and 
subregions with the first and second digits respectively, and follow definitions based on those of Zwally and Giovinetto (201 1 ). The yellow 
border delineates the 2000 m elevation cut-off used in constraining the mascon solution, (b) The AIS. The drainage systems are sequentially 
labeled from 1 to 36. The 2000 m elevation line defining the constraint regions is again delineated by a yellow border. The three Antarctic 
regions of EAIS (drainage systems 3-15, 26-36), WAIS (drainage systems 1, 2, 16-22) and AIS peninsula (drainage systems 23-25) are 
outlined in red. (Note that the yellow elevation line is concealed by the red line along a portion of the boundary between EAIS and WAIS.) 
(c) The GOA mascons. This region is not further delineated by drainage systems or constraint regions. 


Figure 1 shows the GIS, AIS and GOA mascons as well as 
the constraint region boundaries. Also shown in this figure 
are ice-sheet drainage system divides based on those of 
Zwally and Giovinetto (2011). These divides are used for 
computing mass change as the sum of the mascon time 
series within a drainage system, but are not used in the 
mascon estimation. 

We use a characteristic length of 100 km and a 
characteristic time of 10 days in computing the exponential 
constraints for the mascon signal auto-covariance matrix. 
These values represent approximately one parameter 
sample in time and space. We use a damping parameter 
{/j, in Eqn (5)) of 30000 -1 , which provides a similar solution 


contribution from the signal and data covariance to that 
determined from an analysis of the GIS, AIS and GOA 
mascon matrix diagonals. We investigated the impact of the 
damping parameter by performing solutions using values 
that span an order-of-magnitude variation, centered on our 
selected damping parameter value. We find that the 
impact on the estimated mass balance and annual ampli- 
tude is <8% on average over the GIS, AIS and GOA 
mascons, wh i le the 1 a error estimates of the 1 0 day mascon 
solutions can change by as much as 31% on average over 
those regions. Therefore, varying the signal covariance 
damping parameter has a larger impact on the solution 
noise than the signal. 
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Fig. 2. Daily KBRR residual rms solution differences. The daily 
(dots) and 15 day Gaussian smoothed (curves) of the KBRR residual 
rms differences between the v08 solution and the v09 (blue), vIO 
(purple), vl 1 (green) and v12 (orange) solutions. Negative values 
represent a reduction in the residuals and an improvement in the 
forward modeling, including an iteration of the mascon solution. 


Resolution and leakage 

We further analyze the impact of our signal covariance 
(including our choice of the damping parameter and 
characteristic length and time) by performing signal leakage 
analysis within and across constraint regions using the 
resolution matrix identified in Eqn (8). First, we look at the 
influence of the signal covariance on individual mascons. 
Several CIS and AIS mascons were chosen as single-source 
mascons. For each single-source test we fill the mascon 
parameter vector with a time series of mass change for that 
mascon only (we use the time series from our final solution) 
and zero all other mascons in space and time. The resolution 
matrix then multiplies the vector of mascons in time and 
space, where only one mascon contains a time series of 
mass change. Therefore, before application of the resolution 
matrix, only one mascon out of all the global mascons 
contains a nonzero time series of mass change. FHowever, 
after application of the resolution matrix, several neighbor- 
ing mascons will contain nonzero time series, as the signal 
covariance spreads out the signal from the source mascon. 
The effects of the signal covariance are quantified by 
differencing the full mascon vector before and after the 
application of the resolution matrix. The results of several 
tests show that the spatial signal from a single mascon is 
~100% accounted for by adding all the mascons within a 
600 km radius of the original source mascon. Therefore, a 
single mascon has negligible influence at distances greater 
than 600 km. This is equivalent to a Gaussian spatial 
smoothing filter with a 300 km radius. However, the 
constraint region boundaries modify the spatial pattern of 
signal spreading as the signal covariance, by design, limits 
signal leakage across these boundaries. The spatial leakage 
pattern of a single mascon's signal is a function of its 
location with respect to the constraint region boundaries. An 
error analysis of signal leakage across constraint boundaries 
is performed for each region considered. The regional 
leakage analysis and its results are discussed further below. 


In summary, the spatial resolution of our solutions is 
equivalent to a 300 km Gaussian spatial smoothing filter 
within constraint regions, but with significantly reduced 
leakage across constraint region boundaries. 

MASCON SOLUTION ITERATION AND KBRR 
RESIDUAL ANALYSIS 

As described above, we apply a mascon technique to 
estimate global mass change directly from the GRACE KBRR 
observations themselves, taking into account the full Stokes 
noise covariance and a signal covariance. Therefore, we can 
use the reduction in the GRACE KBRR residuals as a 
quantitative objective measure of performance. The KBRR 
residuals are the difference between the GRACE observed 
KBRR data and KBRR data computed from our forward 
models and mascon solution using the NASA Goddard 
Space Flight Center (GSFC) GEODYN precision orbit 
determination and geodetic parameter estimation software. 
Unique to this study, we present the improvement in the 
KBRR residuals from fully iterating the global mascon 
solution (as in Eqn (6)). Because there are >1 1 x 1 0 6 mascon 
parameters, we obtain the solution using iterative methods. 
Following Sabaka and others (2010), we use the precondi- 
tioned conjugate gradient method. Figure 2 presents the 
KBRR residual performance for each solution defined in 
Table 1. This figure shows the difference in the daily root 
mean square (rms) of the KBRR residuals using the v08 
forward model and the KBRR residuals computed from 
subsequent versions of forward models (which include 
iterated mascon solutions for vl 0-v12). Negative differences 
represent improvements due to the new forward model of 
mass change. The v09-v08 daily rms difference in Figure 2 
demonstrates significant reduction in the KBRR residuals 
(model improvement) from forward-modeling hydrology and 
GIA. The residual reduction is seasonal, with the largest 
reduction in daily residual rms in Northern Hemisphere 
spring and an increase in residuals during a portion of 
Northern Hemisphere summer. The residuals show that the 
GLDAS hydrology model improves the global surface mass 
modeling during most of the year, except the Northern 
Hemisphere summer. 

The global mascon solutions represent a reduction in the 
daily KBRR residual rms, as shown in Figure 2. The first 
iteration mascon solution, vIO forward model, represents a 
factor of 2.4 reduction in residuals over the v09 forward 
model (v08 plus GLDAS hydrology and GIA), and further 
iterations represent 5% (vl 1 forward model) and 2% (v12 
forward model) reduction over the preceding iteration. The 
global mascon solutions show reduction in seasonal and 
long-term trend signals in the daily KBRR residuals, due to 
improved modeling of the global mass change, dominated 
by hydrology and land ice. The spatial distribution of the 
residual reduction is presented in Figure 3 as the trend and 
annual amplitude of the GRACE average time-differenced 
KBRR residuals (KBRR-dot). The significant trends in the 
KBRR-dot residuals correspond to the GOA, AIS and GIS 
land-ice regions, and in particular, the largest residual trends 
are found in the West Antarctic ice sheet (WAIS) Amundsen 
and Bellingshausen coasts, and the southeast and northwest 
coasts of the GIS (Fig. 3a). Figure 3a demonstrates the 
significant reduction in these KBRR-dot residual trends using 
the v12 forward model, due to the inclusion of the vl 1 
iterated global mascon solution. 
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Fig. 3. Global maps of the (a) trend and (b) annual amplitude of the time-differenced KBRR residuals (KBRR-dot) for the forward models used 
in the v08 (row 1), v09 (row 2) and v12 (row 3) mascon solutions (Table 1). The differences between v09 and v08 show the impact of the 
additional forward models included in the v09 solution, while the differences between v12 and v09 show the large reduction in residuals 
that is achieved with the iterated mascon solution. 


Iterating the solution has an important impact on the 
regional terrestrial ice mass signals obtained from the global 
mascon solution (Fig. 4; Table 2). Iteration results in 
increases in ice mass loss trend, annual and acceleration 
signal computed from the global mascon solution. It 
accounts for a -16.5, -21.8 and —1.1 Gta -1 correction to 
the GRACE observed CIS, AIS and GOA mascon solution 
trends; and 1 0.7, 34.2 and 4.4 Gt corrections to the GIS, AIS 
and GOA annual amplitude (Table 2). Iteration of the global 
mascon solution is shown to reduce KBRR residuals and 
increase signal recovery. Therefore, we regard the increased 
signal recovery through iteration as an important part of the 
global mascon solution, without which the regional land-ice 
mass signals estimated from the mascon solution would be 
in error by the amounts noted in Table 2 (v12-v09). 

SIGNAL AND ERROR ANALYSIS 

For each of the 41 168 mascons in our global solution, we 
obtain a time series of mass anomalies with 10 day 
sampling. To compute the time series of mass anomalies 
for any region of interest we sum the individual time series at 
each 10 day sample for all mascons contained within the 
region. The 10 day mascon solution ler error is estimated 
from the rms of the misfit of the mascon solution time series 


and a 10 day width Gaussian filter applied to the time series, 
in a similar way to Luthcke and others (2008). Also, as in 
Luthcke and others (2008), the mascon time-series trend and 
annual amplitude are determined from a least-squares 
simultaneous estimate of the trend, annual, semi-annual 
and 1 61 day periodic. The 1 61 day periodic compensates for 
the GRACE alias of an S 2 tide error (Ray and Luthcke, 2006). 
For the trend, we replace the formal error estimate (variance 
of the least-squares parameter) with one that considers a 
first-order autoregression structure for the mascon time 
series (Lee and Lund, 2004). The formal estimates of the 
other parameters (e.g. annual amplitude) are then scaled by 
the ratio of the trend error with autoregressive structure 
considered, and the simple formal trend error. The total error 
is then computed as the root sum square (rss) of the 
following error estimate sources: statistical; leakage in and 
out; atmosphere/ocean forward modeling; and GIA. These 
errors are discussed below. 

Leakage errors 

Leakage of signals in and out of regions of interest occurs 
due to the fundamental spatial resolution limitations of the 
GRACE data themselves and the properties of our mascon 
signal auto-covariance. Errors arising from leakage of signals 
into ('leakage-in') a region of interest (e.g. GIS or AIS) from 
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Fig. 4. The effect of iteration on the mass time series for (a) CIS, (b) AIS and (c) GOA. The third and final iteration solution is shown (v12: 
gray), along with the differences to the original (v09: blue), the first iteration (vIO: purple) and the second iteration (vl 1: green) solutions. 


unmodeled surrounding mass variation sources, such as 
ocean, tides and hydrology, are estimated, as well as land- 
ice mass signal leakage out ('leakage-out') of our regions of 
interest. It should be noted that the hydrology signal leakage 
error has been significantly reduced through the forward 
modeling of a TWS model, as previously discussed, and as 
in Luthcke and others (2008) and Sabaka and others (2010). 
We use a representative global mass anomaly model in 
space and time to estimate the leakage-in and leakage-out 
errors. For this model, we use the vector h/, (Eqn (8)) from 
our final perferred v12 global mascon solution. To estimate 
leakage-in for a region of interest, we zero the region 
mascon parameters in hj. from Eqn (8), but leave the 
remaining global mascon parameter values to produce 
h*/?(o). The estimate of leakage-in error is then the region's 
mascon time series after applying the resolution operator to 
our model, h^( 0 ). We then perform this procedure for each 
region of interest. The leakage-out errors are computed using 
a similar procedure, but in this case all mascon parameters 

Table 2. Effect of solution iteration on mascon regional mass change 
statistics over the time period 1 December 2003 to 1 December 
2010. The reported statistics describe the differences of the final 
and third iteration solution (vl 2) to the original global solution 
(v09), the first iteration solution (vIO) and the second iteration 
solution (vl 1) 


Solution iteration 
by region 

Annual 

amplitude 

Gt 

Acceleration 

Gta^ 2 

Trend 
Gt a -1 

GIS 

vl 2-v09 

10.7 

-0.7 

-16.5 

vl 2-vl 0 

6.3 

-1.0 

-7.8 

v12-v1 1 

2.8 

-0.5 

-3.2 

AIS 

v12-v09 

34.2 

-11.3 

-21.8 

v12-v10 

18.3 

-5.8 

-10.3 

vl 2-vl 1 

7.6 

-2.5 

-4.2 

GOA 

vl 2-v09 

4.4 

2.9 

-1.1 

v12-v10 

1.9 

1.6 

-0.9 

v12-vl 1 

0.8 

0.6 

-0.5 


outside the region of interest are set to zero. The estimate of 
leakage-out error is then the difference between the region's 
model mascon parameters and those after applying the 
resolution operator. A summary of the leakage errors for the 
major land-ice regions of interest is given in Figure 5a and 
Table 3. While the leakage analysis does consider sources 
(1 arcdeg mascons) significantly smaller than the funda- 
mental spatial resolution of GRACE (approximately equiva- 
lent to 300 km radius Gaussian smoother), the analysis did 
not consider point sources and localized signal smaller in 
spatial extent than an individual mascon. An assessment of 
the impact of sub-mascon sources on our leakage estimates 
will need to be evaluated in future analysis. 

Atmosphere/ocean forward-modeling errors 

Following Elan and others (2004), we estimate the errors in 
forward-modeling atmosphere and ocean mass as the 
difference between two sets of models: (1) ECMWF and 
OMCT vs (2) atmospheric mass variations derived from 
NASA/GSFC Modern Era Retrospective-analysis for Research 
and Applications (MERRA) model and a simple ocean 
inverted barometer (IB). A summary of the impact of these 
models on our land-ice solutions is provided in Figure 5b 
and Table 4. The two models are not completely inde- 
pendent, so we do not scale the differences by 1 /\/2 (Elan 
and others, 2004). Ocean mass variation and ocean tides 
forward-modeling error contributions to our land-ice mass 
solutions are included through the leakage-in analysis 
discussed above. The mass anomaly model used in the 
leakage analysis includes the ocean mass and ocean tides 
modeling error observed by GRACE, as the mascon solution 
reflects the mass anomalies from the forward models. 

GIA errors 

The land-ice mass trends estimated from our global mascon 
solutions contain errors due to the uncertainties in the GIA 
model applied. We estimate the GIA error as the difference 
between the minimum and maximum GIA correction over a 
range of plausible GIA models divided by 2 for each region 
of interest. For the GIS and GOA the models considered 
include the nominal ICE-5G Paulson incompressible two- 
layer model discussed previously; an ICE-5G compressible 
model (Geruo and others, 2013); and models based on 
'ANU' ice deglaciation history (Fleming and Lambeck, 
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Fig. 5. (a) The v12 mascon solution mass time series for CIS, AIS and GOA with leakage-in (black circles) and leakage-out (red diamonds) 
estimated errors, (b) The v12 mascon solution time series for GIS, AIS and GOA with atmosphere and ocean model signal (ECMWF/OMCT: 
black circles) and estimated model errors (ECMWF/OMCT - MERRA/IB: red diamonds), (c) The GIS v12 mascon solution time series and the 
difference to the mass time series derived from the averaging kernel filter technique (red diamonds). 


2004) applying a suite of Earth models (27) over a range of 
lithospheric thicknesses, and upper mantle and lower 
mantle viscosities (Barletta and others, 2012, computed 
and provided by V.R. Barletta). For the AIS the models 
considered include the nominal IJ05 R2 used in this study 
and discussed previously; models based on IJ05 R2 with 
varying Earth models; and models based on W12a with the 
'best'-fit Earth model and upper- and lower-bound Earth 
models that give a good fit to the relative sea-level data 
(Whitehouse and others, 2012, computed and provided by 
P. Whitehouse). This is a similar approach to estimate the 
GIA uncertainty to that used in the IMBIE study (Shepherd 
and others, 2012). The post-LIA GIA signal uncertainty is 
±30%, estimated from the model comparisons to GPS 
vertical motions, raised shoreline records of RSL and tide- 
gauge rates of RSL (Larsen and others, 2005). 

Solution technique error analysis 

As another step in our solution evaluation, we compare our 
mascon solution with that of a solution made using the more 
common approach of applying a filter to the monthly 
spherical harmonic fields. In this test we apply a calibrated 


regional GIS averaging kernel to our monthly spherical 
harmonic fields estimated from the same Level-1 B proces- 
sing as our mascon solution (Luthcke and others, 2006b). 
The averaging kernel is constructed, calibrated and applied 
in a similar fashion to Velicogna and Wahr (2005). The 
important aspect of this comparison experiment is that both 
the mascons and the monthly spherical harmonic fields have 
been estimated from the same Level-1 B processing and 
KBRR reduction, so the results isolate the differences arising 
from the technique used to compute the total GIS mass time 
series. Figure 5c shows (in red) the difference between the 
time series computed from the GIS mascons and that 
computed from the averaging kernel; the mascon solution 
time series (black) is also provided in this figure for 
comparison. The differences in the annual amplitude, 
acceleration and trend between the two time series are, 
respectively, 25 Gt, —1 Gta -2 and -11 Gta -1 . These poten- 
tial uncertainties are due solely to the differences in the 
techniques used for the mascon time-series recovery, and 
are on the order of the estimated uncertainties for our 
mascon solution, as summarized in Table 5 and discussed in 
the next subsection. 


Table 3. Effect of leakage-in, leakage-out and the rss of the leakage-in and -out, on the annual amplitude, acceleration and trend of mass 
changes for each main region and subregion over the time period 1 December 2003 to 1 December 2010 


Region 

Annual 

amplitude 

Gt 

Leakage-in 

Acceleration 

Gta- 2 

Trend 

Gta- 1 

Annual 

amplitude 

Gt 

Leakage-out 

Acceleration 

Gta- 2 

Trend 

Gta" 1 

Annual 

amplitude 

Gt 

rss leakage 
Acceleration 

Gta- 2 

Trend 

Gta- 1 

GIS total 

8.7 

-2.9 

2.5 

17.6 

-0.2 

-7.3 

19.7 

2.9 

7.7 

GIS <2000 m 

23.8 

-2.1 

-3.6 

35.8 

1.9 

-17.9 

43.0 

2.8 

18.2 

GIS >2000 m 

15.3 

3.2 

-7.1 

12.2 

1.9 

-2.6 

19.6 

3.7 

7.6 

AIS total 

12.4 

0.5 

-10.9 

16.8 

0.1 

-4.4 

20.8 

0.6 

11.8 

AIS <2000 m 

15.3 

2.9 

-3.6 

16.3 

1.4 

-2.1 

22.4 

3.2 

4.2 

AIS >2000 m 

2.8 

1.0 

-3.3 

5.9 

2.0 

1.7 

6.5 

2.2 

3.7 

AIS east 

6.7 

0.7 

-11.1 

8.6 

1.6 

1.3 

10.9 

1.8 

11.2 

AIS west 

5.1 

-1.7 

-1.7 

5.4 

-3.9 

-8.6 

7.4 

4.3 

8.8 

AIS peninsula 

3.0 

-2.7 

-10.0 

5.3 

-1.8 

-9.1 

6.1 

3.3 

13.5 

GOA total 

19.0 

1.4 

6.9 

34.3 

0.8 

4.2 

39.2 

1.6 

8.1 
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Table 4. Annual amplitude, acceleration and trend of mass changes due to the atmospheric and ocean signal and errors for each main region 
and subregion over the time period 1 December 2003 to 1 December 2010 



Atmospheric 

and ocean signal* 


Atmospheric and ocean error 1 


Region 

Annual amplitude 

Acceleration 

Trend 

Annual amplitude 

Acceleration 

Trend 


Gt 

Gta -2 

Gta' 1 

Gt 

Gta“ 2 

Gta” 1 

GIS total 

199.1 

7.1 

3.5 

9.5 

-2.5 

-4.5 

GIS <2000 m 

90.6 

2.8 

-0.3 

8.9 

-2.0 

-4.1 

GIS >2000 m 

109.9 

4.3 

3.7 

1.7 

-0.5 

-0.4 

AIS total 

697.7 

-58.5 

-5.8 

17.6 

-21.5 

8.0 

AIS <2000 m 

262.9 

-40.4 

-16.9 

7.7 

-13.5 

-6.5 

AIS >2000 m 

437.0 

-18.1 

11.1 

10.3 

-7.9 

14.5 

AIS east 

561.2 

-37.4 

2.2 

15.1 

-17.6 

8.3 

AIS west 

127.6 

-17.7 

-5.9 

8.3 

-2.2 

-0.1 

AIS peninsula 

9.9 

-3.4 

-2.1 

4.2 

-1.7 

-0.2 

GOA total 

42.8 

-0.7 

-2.1 

1.7 

-0.6 

-1.8 


*Sum of ECMWF atmospheric model and OMCT ocean model. 

' I titterence between the combined ECMWF/OMCT models and the combined MERRA/IB models. 


Computation of seasonal balances 

We use a linear regression analysis to extract the multi-year 
average trend, acceleration and annual amplitude from our 
10 day mascon time series. One of the strengths of GRACE 
solutions, and of the mascon solution presented here in 
particular, is the high temporal resolution at which the 
evolution of land-ice mass balance can be observed and 
quantified. Therefore, we seek to compute annual mass 
balances from our GRACE mascon time series, in order to 
better quantify the more complex spatial and temporal 
variability of land ice. Luthcke and others (2008) simply 
used the extrema of the full mascon time series as the start 
and end of the balance years from which to compute the 
annual mass balances. Because the mass time series for the 
sum of GOA glacier mascons has a very clean annual signal 
with little additional signal or noise, this is an acceptable 
approach. However, for the more complex mascon time 


series from the GIS, WAIS, East Antarctic ice sheet (EAIS) and 
Antarctic ice sheet peninsula (AISP), an adaptive filter 
approach is necessary. Therefore, we employ an ensemble 
empirical mode decomposition (EEMD) adaptive filter to first 
compute the timing of the balance seasons (accumulation 
and loss), and then to compute the annual mass balances. 

The EEMD is a noise-assisted data analysis improvement 
of the EMD method, which provides an adaptive decom- 
position of the time series into components known as 
intrinsic mode functions (IMFs) (Wu and Huang, 2009). The 
IMFs are monochromatic at a given point in time, but vary in 
time with both amplitude and frequency. The IMF adaptive 
decomposition separates out the frequency content of the 
time series, which represent different physical modes. An 
ensemble (many runs) of the EMD is performed, applying a 
new generation of noise to the original signal on each run 
within the ensemble. The added finite white noise provides a 


Table 5. Summary of the vl 2 global mascon solution for the land-ice regions studied over the time period 1 December 2003 to 1 December 
2010. The reported parameters are computed after applying the GIA and LIA model correction to the mascon time series. The total 
uncertainties include contributions from statistical errors, mismodeling of the atmosphere, leakage errors, GIA model errors and, for Gulf of 
Alaska, the LIA model errors 


Region 

Area 
1 0 6 km 2 

Annual amplitude 
Gt 

GRACE mascon solution 
Acceleration 
Gta~ 2 

Trend 

Gta- 1 

GIA model* 
Trend 
Gta- 1 

GIS total 

2.46 

131 .9 ±28.6 

-10.1 ±6.3 

-230.3 ±12.1 

4.7 ± 6.7 

GIS <2000 m 

1.49 

109.0 ±48.6 

-1 0.6± 6.7 

-223. 7± 19.8 

6.9±3.8 

GIS >2000 m 

0.98 

23.1 ±22.3 

0.5 ±4.7 

-6.6 ±8.6 

-2.2 ±3.0 

AIS total 

13.96 

31 6.4 ±57.7 

-32.6±25.6 

-80.8 ±26.2 

47.0 ±17.9 

AIS <2000 m 

6.93 

172.6 ±50.7 

-32.6 ±18.6 

-1 1 3.3 ±21.0 

42. 0± 15.9 

AIS >2000 m 

7.02 

1 43.7 ± 25.1 

0.1 ±10.2 

32.5 ±21.6 

5.0± 14.6 

AIS east 

10.57 

215.4±48.8 

22.3 ±21 .6 

62.8±28.1 

1 6.8 ±21 .6 

AIS west 

2.76 

83.2 ±17.6 

-45.5 ±6.1 

-1 06. 0± 15.9 

26.2 ± 12.8 

AIS peninsula 

0.63 

1 8.3 ±14.6 

-9.3 ±5.1 

-37.7 ± 1 4.0 

4.1 ±1.6 

GOA total 1 

0.76 

165.9±46.2 

1 .7 ± 6.9 

-68.8 ±11.0 

2.1 ±2.0 


*GIS and GOA are corrected with the ICE5G_Paulson GIA model. AIS is corrected with the IJ05_R2 GIA model. 
tThe trend correction for the GOA region from the Larsen LIA model is 10.2±3.1 GtaW 
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means to sort out signals with similar scale into their proper 
IMFs and thus mitigate mode mixing. The mean of the 
ensemble is treated as the final true result, so the effect of 
noise is minimized. Unlike typical signal analysis tech- 
niques, such as the linear regression analysis applied above, 
the EEMD is an adaptive filter and does not rely on linear 
and stationarity assumptions. 

We apply the EEMD to adaptively find the extrema of the 
non-stationary annual signal within the mascon time series of 
interest. The minima and maxima of the annual IMF represent 
the timing of the mass-balance seasons (beginning of the 
mass accumulation and mass loss seasons). This is the most 
important aspect of our EEMD application. The adaptive filter 
is used to isolate the annual signal, and then to identify the 
timing of the balance seasons and length of each balance 
year as the annual signal is non-stationary and evolving. The 
original signal is then reconstructed by summing the IMFs, 
and eliminating noise by discarding the first few IMFs (Fig. 8). 
We then compute the mass value of the reconstructed signal 
at the beginning of the mass-balance seasons (annual IMF 
minima). We use 100 runs within each EEMD ensemble, 
where the added white noise is recomputed on each run 
within the ensemble, with a noise ratio equal to twice the 
ratio of the standard deviation of the mascon solution noise 
(difference between 10 day solution time series values and a 
30 day Gaussian filter), and the standard deviation of the 
signal (with trend removed). The EEMD is then run 50 times, 
where the overall signal trend is removed before each run 
and is added back after the reconstruction of the signal from 
the selected IMFs. The mean (computed over the 50 EEMD 
runs) timing value of the annual IMF minima and corres- 
ponding reconstructed signal mass values are used to 
compute the annual balance (Fig. 8). The annual balance is 
computed by differencing the reconstructed mass values at 
the times of two consecutive mean annual IMF minima. The 
mass difference, or annual balance, is labeled with the year 
that is mostly encompassed between the two minima. 
Therefore, for the GIS and GOA land ice, the balance year 
starts in the late summer to mid-fall of the previous year, 
while for the AIS the balance year starts in the Northern 
Hemisphere spring. The standard deviation of these values 
over the 50 runs is used as an estimate of the noise 
contribution to the uncertainty. To estimate a potential 
systematic error in our balance year timing and values, the 
EEMD analysis results are differenced from an iterative 
piecewise trend and annual linear regression analysis of the 
mascon time series, where the annual pieces are constrained 
to form a continuous signal. The error estimate of balance 
year extrema timing and mass value is then computed as the 
rss of the noise and systematic error estimates. The error 
estimate of the annual balance is then the rss of the error 
estimates for the two balance year extrema. 


LAND-ICE MASS EVOLUTION RESULTS 

In this section we provide a summary of our iterated global 
mascon solution, focusing on the results for the five most 
significant land-ice regions: GIS, WAIS, EAIS, AISP and 
GOA. Table 5 provides an overall summary of the regional 
mass time series computed from the mascon solution. While 
the time series cover a longer time period, the model 
parameter fit is applied to an integer number of years, 
1 December 2003 to 1 December 2010. We use the 
terminology defined by Cogley and others (2011) to explain 



Year 


Fig. 6. The v12 10 day mascon solution mass time series (circles) 
and 1 0 day Gaussian smoothed solution (curves) for GIS total (blue), 
GIS <2000 m elevation (purple) and GIS >2000 m elevation (green). 
The red lines show the best-fit trend plus acceleration. 

glacier and ice-sheet variations in terms of variations in their 
climatic balances, £> c | im , and calving fluxes, D, ignoring the 
contributions from basal balances. 

GIS results 

The time series for the sum of the GIS mascons (GIS total), 
and the mascons above and below 2000 m are shown in 
Figure 6, and a summary of these time series is provided in 
Table 5. The spatial distribution of the GIS trend and 
acceleration is shown in Figure 7. The GIS exhibits the 
largest rates of mass loss relative to both Antarctica and 
Alaska, with a trend of —230 ±12 Gta -1 , an acceleration of 
mass loss of —10 ±6 Gta -2 and an average annual ampli- 
tude of 132±29Gt. Both the magnitude of the GIA 
correction and our estimated uncertainty determined from 
the model differences are quite small in comparison with the 
total corrected trend (Table 5). Our estimates agree to within 
stated error bars with previous estimates spanning at least 
the period 2003 to October 2009: -219±38Gta -1 for 
April 2002-November 2009 (Chen and others, 2011); 
-238 ± 29 Gta -1 for October 2003-October 2009 (Sasgen 
and others, 2012); — 222 ±9 Gta -1 for January 2003- 
December 2010 (Jacob and others, 2012). Our estimated 
acceleration of — 10±6Gta -2 is in agreement with 
-15±7Gta -2 for August 2001 -August 2010 (Sasgen and 
others, 201 2) and — 1 7 ± 8 Gta -2 for April 2002-June 201 0 
(Rignot and others, 2011). The GIS mass evolution is 
dominated by the signal from the coastal margins of the 
ice sheet, as defined in our solution as mascons <2000 m. 
Large coastal mass loss rates result from high rates of D 
associated with the presence of fast-flowing marine-termi- 
nating glaciers along the northwest, southwest and southeast 
regions of the GIS (Moon and others, 2012), and due to the 
higher surface melt rates that occur at low elevations on the 
ice sheet. Above 2000 m there is little statistically significant 
signal. 

The GIS annual balance results from the application of the 
EEMD filter are shown in Figures 8 and 9. The GIS annual 
balance anomalies (difference from the mean) show large 
interannual variability (Fig. 9). The GIS mean annual balance 
over the balance years 2004-10 is — 248 ± 12 Gta -1 . It is 
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Fig. 7. The G1S v12 mascon solution mass changes computed over 
the time period 1 December 2003 to 1 December 2010. (a) Trend 
corrected with ICE5G_Paulson GIA model and (b) acceleration. 


important to note that our mean annual balance is slightly 
different from the trend reported above, primarily due to the 
differing time periods, where integer balance years are 
needed for the annual balance computation. The balance 
years 2005-09 have anomaly magnitudes <20% of the mean 
annual balance and alternate between positive and negative 
values. In contrast, the GIS balance years 2004 and 2010 
have large positive and negative anomalies, ~50% of the 
mean annual balance, giving rise to the statistically signifi- 
cant acceleration determined from the linear regression 
analysis (Table 5). Therefore our reported GIS mass loss 
acceleration is, in part, an artifact of the GRACE sampling 
interval, rather than an indication of a persistent multi-year 
signal, and longer time series will be needed before making 
inferences about future evolution of the ice sheet. 

In the northwest GIS (drainage systems 81 and 82) we 
observe an increase in mass losses after 2005 (Figs 10 and 
1 1), attributed to a reduction in snow accumulation (Sasgen 
and others, 2012) that combines with increased D after 
2007 (Moon and others, 2012) to further enhance mass 
losses. The timing and extent of the northward progression of 
mass loss is confirmed by GPS observations, showing 
enhanced rates of uplift late in 2005 at Thule Station in 


northwest Greenland (Khan and others, 2010). In the 
southwest GIS (drainage systems 61, 62 and 72) we observe 
persistent mass losses, attributed to high D that commenced 
in the late 1990s, that are dominated by the large calving 
losses from Jakobshavn Isbras, the largest marine-terminating 
outlet on the GIS (Joughin and others, 2012). 

Large interannual mass variability is observed in drainage 
system 32, the region of the Geikie Peninsula in central East 
Greenland that is largely composed of ice not receiving flow 
contributions from the GIS (Fig. 10) (Jiskoot and others, 
2012). The mass time series for this region shows an increase 
in mass loss beginning in 2005, followed by mass gains from 
2007 to 09 (Fig. 1 1 ). The glaciers in this area are smaller and 
steeper than ice-sheet outlet glaciers and exist in a region 
with high precipitation rates (Ettema and others, 2009), so 
they likely respond rapidly to climate variations in a fashion 
similar to other maritime glacier systems. Approximately 
90% of the glaciers in this region drain through tidewater 
glacier outlets (Jiskoot and others, 2012), and those draining 
to the south into the Blosseville Kyst may be particularly 
sensitive to changes in climate, due to their southerly aspect 
and concentration of ice at lower elevations. 

Several studies have reported synchronous retreat of 
marine-terminating glaciers in southeast Greenland during 
2000-06 (Howat and others, 2008), followed by advance 
and thickening of Helheim and Kangerdlugssuaq Glaciers 
(Howat and others, 2007; Joughin and others, 2008). The 
reduction in mass loss associated with these advances has 
been assumed to be synchronous across the southeast coast 
(Murray and others, 2010), although evidence for this was 
based, in part, on terminus advance rates that did not 
measure glacier mass or volume changes directly (Howat 
and others, 2008). Our results show a deceleration in mass 
loss during the entire 2003-10 period (Fig. 7b), but analysis 
of the individual mass-balance seasons shows that much of 
the deceleration occurred due to mass gains in the central 
east GIS in 2008, and a synchronous reduction in mass 
losses in 2009 (Fig. 10). These findings are consistent with 
those of Chen and others (2011), who used GRACE 
observations to infer a significant reduction in the mass loss 
of southeast GIS glaciers during September 2007-November 
2009; however, our results pinpoint the timing of this 
reduction in mass loss to the 2009 balance year. The 2009 
slowdown of mass loss in the southeast GIS likely occurred 
due to a complex array of factors, including a deceleration 
and advance of several southeast GIS glaciers (Murray and 
others, 2010), as well as 2008 and 2009 positive accumu- 
lation anomalies (Sasgen and others, 2012). We note that 
observed deceleration and advance of southeast GIS glaciers 
occurred 1-2 years prior to our observed 2009 reduction in 
mass losses. This reflects the complex and asynchronous 
nature of ice-sheet dynamics in the southeast GIS (Moon and 
others, 2012). 

Our observations show evidence of the response of the 
GIS to variations in summer air temperatures and the length 
of the summer melt season. We observe high rates of mass 
loss across the western and northwestern GIS during 2008, 
associated with extreme summer snowmelt that summer, 
detected by passive microwave imaging sensors (Tedesco 
and others, 2008). Coastal mass balances during 2009 were 
less negative than all other years in the GRACE record, and 
this was followed in 2010 with the most negative mass 
balances below 2000 m (Fig. 10). These patterns agree well 
with observed and simulated melt extent records showing 
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Fig. 8. The EEMD v12 mascon solution time series and error analysis for (a) CIS, (b) GOA, (c) EAIS, (d) WAIS and (e) AIS peninsula. Each plot 
shows the v12 mascon solution time series (black circles), the mean of 50 EEMD runs, where each case sums the IMFs from the one 
preceding the 'annual' IMF to the final IMF (blue curve), the mean values from the 50 EEMD runs at the times of the 'annual' IMF extrema 
(red triangles) and the extrema error ellipses (green ellipses). 


low/high melt extents in the 2009/1 0 summers (Mernild and 
others, 201 1). 

AIS results 

The mascon solution for the AIS regions is summarized in 
Table 5 and Figures 12 and 13. Unlike the CIS, both the 
magnitude and uncertainties of the GIA corrections for the 
AIS are significant with respect to the total trend. The overall 
AIS trend for the December 2003-December 2010 time 
period studied is — 81±26Gta -1 . Our trend estimate 
agrees, within stated uncertainties, with -69±18Gta -1 
obtained by King and others (2012) for August 2002- 
December 2010, using a comparable recent regional GIA 
model. The AIS trend is dominated by significant mass loss 
from the WAIS, with a trend of —1 06 ± 1 6 Gta -1 . A trend of 
— 38 ± 14 Gta -1 is found for the AIS peninsula, while the 
EAIS trend is 63 ± 28 Gta -1 . The largest trends are found in 
the WAIS along the Amundsen and Bellingshausen Sea 
coasts (drainage systems 19-22), concentrated at the Pine 


Island embayment (convergence of drainage systems 
19-22), and at the northern tip of the peninsula (drainage 
system 24). These results are in general agreement with 
— 132 ± 60 Gta -1 for the WAIS and —60 ± 46 Gta -1 for the 
AIS peninsula, determined from radar interferometry and 
regional climate models for the year 2006 (Rignot and 
others, 2008). In addition, our AIS peninsula trend is in good 
agreement with — 42 ±9 Gta -1 for January 2003-March 
2009, determined from a GRACE spherical cap mascon 
solution (Ivins and others, 2011). 

The most striking signal, with the greatest potential 
consequence for sea-level rise, is the large accelerated mass 
loss exhibited by the WAIS at —46 ± 6 Gta -2 . While we do 
observe interannual variability of the WAIS annual balances, 
there is a consistent negative trend in these data, repre- 
senting a persistent multi-year acceleration of mass loss 
(Figs 8d and 9). The largest accelerated loss is concentrated 
at Pine Island embayment and along the Amundsen and 
Bellingshausen Sea coasts with significant mass loss 
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Fig. 9. The annual mass balance mean (2004-10) and annual anomalies about the mean for GIS (blue), GOA (green), EAIS (purple), WAIS 
(orange) and AIS peninsula (cyan) from the EEMD analysis of the v12 mascon solution. The values and corresponding error bars are 
determined from the mean extrema and extrema error ellipses shown in Figure 8. 
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Fig. 10. The CIS annual mass balances and mean annual mass balance determined from the EEMD analysis of the vl 2 mascon solution. The 
timing of the seasons is determined by the mean EEMD extrema in Figure 8a, and the mass values are determined by the result of the EEMD 
method (summing the IMFs from the one preceding the 'annual' IMF to the final IMF) for each individual mascon time series at the time of 
the total GIS extrema. 


extending to the Zumberge coast (drainage system 1) during 
the 2009 balance year (Fig. 14). These changes are consist- 
ent with observations of grounding-line retreat (Rignot and 
others, 2008), acceleration (Joughin and others, 2010) and 
drawdown (Wingham and others, 2009) of the Pine Island 
Glacier driven by basal melting of its large ice shelf 
(Pritchard and others, 2012). 

The EAIS trend (63±28Gta _1 ) and accelerated mass 
gain (22 ±22 Gta -2 ) (Table 5; Fig. 8c) are not persistent 
multi-year signals but are due to an exceptional positive 
surface mass-balance anomaly during the 2009 accumu- 
lation season, concentrated primarily along the coast of 
Dronning Maud Land (drainage systems 5-7), but also along 
the coast of Wilkes Land (drainage systems 11 and 12) 
(Figs 8c and 14). These observations have similar patterns to 
melting days anomalies derived from passive microwave 
observations which, for the entire AIS, reached a new 
historical minimum for the 1980-2009 period (Tedesco and 
Monaghan, 2009). In addition, our GRACE observations 
agree quite well with the timing and spatial location of a 
significant 2009 accumulation anomaly observed in the 
surface mass-balance record derived from the regional 
atmospheric climate model RACM02 (Lenaerts and others, 
2012; Shepherd and others, 2012). 

The northern tip of the AIS peninsula exhibits consistent 
negative annual balance, with the exception of the 2008 
balance year (Fig. 14), which gives rise to the shallow 


positive acceleration shown in Figure 13. The southern area 
of the peninsula shows small positive and negative annual 
balances, with the exception of large negative annual 
balances for 2007 and 2009, giving rise to the sharp losses 
shown in the AIS peninsula time series for these years 
(Fig. 8e). These large negative annual balance anomalies 
give rise to the acceleration signal in the southern peninsula 
shown in Figure 13, and the overall peninsula acceleration 
of — 9±5Gta -2 . Large mass losses on the AIS peninsula 
result from tributary glaciers that feed into ice shelves that 
have recently disintegrated (Shuman and others, 201 1), and 
increasingly negative mass balances resulting from strong 
climatic warming in the region. Our observations agree with 
satellite observations of surface lowering (Fricker and 
Padman, 2012) and model estimates of increasing mass 
losses from peripheral glaciers and ice caps in the region 
(Flock and others, 2009). 

GOA results 

The mascon solution for the GOA glacier region is 
summarized in Table 5 and Figures 8b and 15. The overall 
trend of the GOA glaciers for the time period studied is 
—69 ± 1 1 Gta -1 , and —77 ± 1 1 Gta -1 if we take the mean 
of the 2004-10 balance years (Fig. 9). Most of the GIA 
correction in this region results from uplift corrections due to 
LIA mass loss in the Glacier Bay region of Alaska. Our 
estimates very nearly agree to within stated error bars with 
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Fig. 11. The v12 10 day mascon solution time series (circles) and 
10 day Gaussian smoothed solution (curves) for selected drainage 
systems within the CIS and AIS, and mascons within the GOA land- 
ice region. The colors of the plotted curves correspond to the 
drainage system and mascon colors in Figure 1. 


-46±7Gta^’ reported by Jacob and others (2012). 
Differences in the best estimates between these two studies 
likely results from different treatment of TWS, and slightly 
different solution domains. 

GOA glaciers exhibited large seasonal and interannual 
variability during the 2003-10 measurement period (Figs 8b, 
11 and 16). We observed large mass losses during summer 
2004, in response to record high temperatures across Alaska 
that year (Truffer and others, 2005). Even larger mass losses 
occurred in 2009, despite the absence of a strong tempera- 
ture forcing, and we attribute these losses to the 31 March 
2009 eruption of Mount Redoubt that spread volcanic ash 
across much of the GOA region (Schaefer and others, 2012) 
and likely enhanced melt rates through a reduction in surface 
albedo. Large 2007/08 winter accumulation was followed by 
a cool 2008 summer season, resulting in near-balance 
conditions during 2008. This offset the strong 2004 and 
2009 melt seasons and resulted in no significant acceleration 
in GOA mass balance during the GRACE period. 

The spatial distribution of mass balances largely follows 
the distribution of glacierization in the GOA region, with 
largest mass losses occurring in the St Elias and Glacier Bay 
regions. The highly variable distribution of glaciers in the 
GOA region, with large signals occurring in the center of the 
domain and much smaller signals at the northwestern and 
southeastern edges, likely results in smearing and leakage of 
signal between individual mascons. Therefore, although the 
total GOA signal is not affected by these problems, 
additional processing steps are required to further constrain 
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Fig. 12. The v12 10 day mascon solution time series (circles) and 
1 0 day Gaussian smoothed solution (curves) for the AIS total (blue), 
AIS <2000 m elevation (purple), AIS >2000 m elevation (green), 
EAIS (orange), WAIS (cyan) and AIS peninsula (gold). The red curves 
show the best-fit trend and acceleration. 



Fig. 13. The AIS v12 mascon solution mass changes computed over 
the time period 1 December 2003 to 1 December 2010. (a) Trend 
corrected with the IJ05_R2 GIA model and (b) acceleration. The 
EAIS, WAIS and AIS peninsula regions are demarcated by the cyan 
borders. 
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Fig. 14. The AIS annual mass balances and mean annual mass balance determined from the EEMD analysis of the vl 2 mascon solution. The 
timing of the seasons is determined by the mean EEMD extrema in Figure 8c, d and e, and the mass values are determined by the result of the 
EEMD method (summing the IMFs from the one preceding the 'annual' IMF to the final IMF) for each individual mascon time series at the 
time of the extrema for the region in which the mascon is located (east, west or peninsula). The EAIS, WAIS and AIS peninsula regions are 
demarcated by the cyan borders. The mean annual mass balance is shown for the v12 solution corrected with the IJ05_R2 CIA model. 




Fig. 15. The GOA mascon trends computed from the v12 mascon 
solution corrected with the ICE5G_Paulson GIA model and the 
Larsen LIA model over the time period 1 December 2003 to 
1 December 2010. 


mass changes over drainage systems smaller than our 
current solution domain. 

CONCLUDING REMARKS 

We have developed a new global mascon solution spe- 
cifically tuned to observe high-latitude land-ice mass 
evolution estimated directly from the GRACE inter-satellite 
range-rate data. The study is unique in many ways including: 
(1) a detailed error analysis of a constrained global mascon 
solution with estimates of signal leakage in and out of the 
regions of interest; (2) a controlled comparison to solutions 
derived from monthly spherical harmonic solutions using 
the same GRACE data processing and forward models, 
providing estimates of uncertainty due to solution technique; 
(3) fully iterating a global mascon solution to convergence 
and quantifying the signal recovery improvement and 
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Fig. 16. The GOA land-ice region annual mass balances and mean annual mass balance determined from the EEMD analysis of the v12 
mascon solution. The timing of the seasons is determined by the mean EEMD extrema in Figure 8b, and the mass values are determined by 
the result of the EEMD method (summing the IMFs from the one preceding the 'annual' IMF to the final IMF) for each individual mascon time 
series at the time of the total GOA extrema. 
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GRACE inter-satellite range-rate residual reduction; and 
(4) applying the EEMD to compute land-ice annual 
balances. We estimate the land-ice mass trends for the time 
period 1 December 2003 to 1 December 2010 to be 
-230 ± 12 Gta -1 for the GIS, -81 ± 26 Gta -1 for the AIS 
and — 69±11Gta -1 for the GOA glaciers. For the 
subregions of the AIS we estimate 63 ±28 Gta -1 for the 
EAIS, — 1 06 ±16 Gta -1 for the WAIS and -38 ± 14 Gta" 1 
for the AIS peninsula. The total over the land-ice regions 
studied is -380 ± 31 Gta -1 , equivalent to -1.05 ±0.09 
mm a -1 sea-level rise. Over the same time period we 
estimate the mass acceleration to be -1 0 ± 6 Gta -2 for the 
GIS, —33 ±26 Gta -2 for the AIS and 2±7Gta -2 for the 
GOA glaciers, giving a total of -41 ± 27 Gta -2 , equivalent 
to —0.1 1 ± 0.08 mm a -2 sea-level rise. For the subregions of 
the AIS we estimate the accelerations to be 22 ± 22 Gta -2 
for the EAIS, —46 ± 6 Gta -2 for the WAIS and — 9 ± 5 Gta -2 
for the AIS peninsula. The trends and accelerations deter- 
mined over our study period have been shown to be highly 
dependent on significant seasonal and annual balance 
anomalies, making it difficult to predict future land-ice 
mass balance and its contribution to sea level. One possible 
exception is the WAIS acceleration signal, which shows a 
persistent negative trend in annual balance. The iterated 
global mascon solution presented here provides the spatial 
and temporal monitoring of land-ice evolution to aid in 
improved understanding and modeling. 
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